{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Lenovo\\AppData\\Roaming\\Python\\Python311\\site-packages\\tqdm\\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from mgwr.sel_bw import Sel_BW\n",
    "from mgwr.gwr import MGWR\n",
    "import pandas as pd\n",
    "import shap\n",
    "import numpy as np\n",
    "import shap\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 读取CSV文件\n",
    "data = pd.read_csv(r\"D:\\Lenovo\\Desktop\\云南大学\\毕业设计\\毕设数据\\输出数据\\csv\\countyzhujiang_3\\countyzhujiang_Factors_2020_3.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 提取坐标\n",
    "coords = np.array(list(zip(data['X_position'], data['Y_position'])))\n",
    "\n",
    "# 因变量y\n",
    "y = data['总碳排2020'].values.reshape((-1, 1))\n",
    "\n",
    "# 自变量X\n",
    "X_original = data[['SMCI_2020', 'GDP_2020', '人口_2020','土地利用结构_2020'\n",
    "          ,'人口密度_2020','年均温_2020','人均GDP_2020','年降水量_2020',\n",
    "          '燃烧面积占比_2020_2','火点个数_2020','单位面积FRP_2020']].values\n",
    "scaler=StandardScaler()\n",
    "X=scaler.fit_transform(X_original)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 设置每个自变量的最小/最大带宽（以相邻要素数为单位）\n",
    "# 注意每年的带宽搜索范围不同，带宽范围根据“指标体系”表±10得到\n",
    "# 截距项带宽也包含在其中，截距是第一个x项，截距项带宽搜索范围为全部\n",
    "bw_min = [10,10, 20, 25,20,16,10,20,10,60,16,50]\n",
    "bw_max = [350,350, 40, 45,40,36,350,40,350,80,36,70]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Backfitting: 100%|██████████| 100/100 [12:53<00:00,  7.73s/it]\n"
     ]
    }
   ],
   "source": [
    "# 初始化带宽选择器（启用多尺度模式）\n",
    "# Sel_BW是带宽类\n",
    "selector = Sel_BW(coords, y, X, multi=True, fixed=False,n_jobs=-1,\n",
    "                  constant=True)\n",
    "\n",
    "# 搜索最优带宽（使用AICc准则）\n",
    "bw=selector.search(\n",
    "    verbose=False,\n",
    "    criterion='AICc',\n",
    "    multi_bw_min=bw_min ,\n",
    "    multi_bw_max=bw_max,\n",
    "    max_iter=100,\n",
    "    max_iter_multi=100,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Inference: 100%|██████████| 16/16 [00:00<?, ?it/s]\n"
     ]
    }
   ],
   "source": [
    "# 初始化模型（使用自适应带宽）\n",
    "mgwr_model = MGWR(\n",
    "    coords, y, X,\n",
    "    selector=selector,\n",
    "    fixed=False,  # 使用k邻近（自适应带宽）\n",
    "    kernel='gaussian',\n",
    "    hat_matrix=True,\n",
    "    n_jobs=-1,\n",
    ")\n",
    "\n",
    "# 拟合模型\n",
    "results = mgwr_model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===========================================================================\n",
      "Model type                                                         Gaussian\n",
      "Number of observations:                                                 357\n",
      "Number of covariates:                                                    12\n",
      "\n",
      "Global Regression Results\n",
      "---------------------------------------------------------------------------\n",
      "Residual sum of squares:                                       127064950339123249152.000\n",
      "Log-likelihood:                                                   -7720.370\n",
      "AIC:                                                              15464.740\n",
      "AICc:                                                             15467.801\n",
      "BIC:                                                           127064950339123249152.000\n",
      "R2:                                                                   0.864\n",
      "Adj. R2:                                                              0.860\n",
      "\n",
      "Variable                              Est.         SE  t(Est/SE)    p-value\n",
      "------------------------------- ---------- ---------- ---------- ----------\n",
      "X0                              814777990.573 32119533.568     25.367      0.000\n",
      "X1                              173235095.224 58585331.633      2.957      0.003\n",
      "X2                              -297500044.113 69191225.743     -4.300      0.000\n",
      "X3                              1666616287.502 64467276.230     25.852      0.000\n",
      "X4                              386417374.211 77890443.858      4.961      0.000\n",
      "X5                              -477419767.825 67702489.207     -7.052      0.000\n",
      "X6                              -64901608.722 58597728.279     -1.108      0.268\n",
      "X7                              322995232.428 52657693.929      6.134      0.000\n",
      "X8                              -36838009.155 41090835.109     -0.897      0.370\n",
      "X9                              63031208.636 42178509.236      1.494      0.135\n",
      "X10                             -199320358.609 47473462.366     -4.199      0.000\n",
      "X11                             88328443.458 47045607.852      1.878      0.060\n",
      "\n",
      "Multi-Scale Geographically Weighted Regression (MGWR) Results\n",
      "---------------------------------------------------------------------------\n",
      "Spatial kernel:                                           Adaptive gaussian\n",
      "Criterion for optimal bandwidth:                                       AICc\n",
      "Score of Change (SOC) type:                                     Smoothing f\n",
      "Termination criterion for MGWR:                                       1e-05\n",
      "\n",
      "MGWR bandwidths\n",
      "---------------------------------------------------------------------------\n",
      "Variable             Bandwidth      ENP_j   Adj t-val(95%)   Adj alpha(95%)\n",
      "X0                     153.000      1.327            2.086            0.038\n",
      "X1                     348.000      1.058            1.991            0.047\n",
      "X2                      39.000      1.840            2.218            0.027\n",
      "X3                      26.000      4.629            2.563            0.011\n",
      "X4                      39.000      2.695            2.365            0.019\n",
      "X5                      35.000      2.366            2.316            0.021\n",
      "X6                     349.000      1.070            1.996            0.047\n",
      "X7                      39.000      3.429            2.455            0.015\n",
      "X8                     349.000      1.098            2.007            0.046\n",
      "X9                      79.000      2.489            2.335            0.020\n",
      "X10                     35.000      4.137            2.523            0.012\n",
      "X11                     69.000      2.252            2.297            0.022\n",
      "\n",
      "Diagnostic information\n",
      "---------------------------------------------------------------------------\n",
      "Residual sum of squares:                                       40239945064092418048.000\n",
      "Effective number of parameters (trace(S)):                           28.390\n",
      "Degree of freedom (n - trace(S)):                                   328.610\n",
      "Sigma estimate:                                                349935943.073\n",
      "Log-likelihood:                                                   -7515.124\n",
      "AIC:                                                              15089.028\n",
      "AICc:                                                             15094.498\n",
      "BIC:                                                              15202.996\n",
      "R2                                                                    0.957\n",
      "Adjusted R2                                                           0.953\n",
      "\n",
      "Summary Statistics For MGWR Parameter Estimates\n",
      "---------------------------------------------------------------------------\n",
      "Variable                   Mean        STD        Min     Median        Max\n",
      "-------------------- ---------- ---------- ---------- ---------- ----------\n",
      "X0                   1018839964.053 246272601.037 702674681.788 1071510276.832 1343604665.766\n",
      "X1                   91425574.844 8581708.645 77785040.281 89653862.395 106804639.437\n",
      "X2                   949077016.155 1039156751.667 -655355261.575 1090026891.026 3147351662.502\n",
      "X3                   1341014915.109 710824065.157 -437635999.831 1351164113.076 2568413138.994\n",
      "X4                   683892792.693 295662933.249 234309711.529 729144853.128 1461176649.370\n",
      "X5                   -935454793.137 684696918.935 -2390620663.993 -671681017.924 108013150.322\n",
      "X6                   -101626831.836 8993951.447 -115969186.735 -102832192.387 -75603701.480\n",
      "X7                   323222713.167 214920327.766 -212054019.042 387738600.376 748016729.612\n",
      "X8                   -6644087.109 9061722.727 -24496743.306 -2474845.327 6286262.702\n",
      "X9                   51275525.565 130733387.402 -177544305.770 79778059.007 264041909.781\n",
      "X10                  -101935441.910 349633666.663 -580067528.081 -209298587.074 899448030.099\n",
      "X11                  53245367.054 83176279.347 -156857638.194 84089237.521 158936014.743\n",
      "===========================================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取变量帽子矩阵S（n*n*k）\n",
    "R=results.R\n",
    "# 确保y和predy为一维数组\n",
    "y_flat = results.y.flatten()\n",
    "predy_flat = results.predy.flatten()\n",
    "\n",
    "# 计算每个点的局部R2\n",
    "local_R2 = []\n",
    "for i in range(R.shape[0]):\n",
    "    # 初始化有效邻域点集合\n",
    "    neighbors = set()\n",
    "    \n",
    "    # 遍历所有变量，合并各变量的局部邻域\n",
    "    for j in range(R.shape[2]):\n",
    "        # 获取第j个变量的权重行\n",
    "        weights_ij = R[i, :, j]\n",
    "        # 找到该变量的非零权重点（对应其带宽内的邻域）\n",
    "        var_neighbors = np.where(weights_ij != 0)[0]\n",
    "        neighbors.update(var_neighbors)\n",
    "    \n",
    "    # 转换为列表\n",
    "    neighbors = list(neighbors)\n",
    "    \n",
    "    if not neighbors:\n",
    "        # 若无邻域点，R²设为0\n",
    "        local_R2.append(0.0)\n",
    "        continue\n",
    "     # 提取邻域内的y和预测值\n",
    "    y_neigh = y_flat[neighbors]\n",
    "    y_hat_neigh = predy_flat[neighbors]\n",
    "    \n",
    "    # 计算SS_res和SS_tot\n",
    "    ss_res = np.sum((y_neigh - y_hat_neigh) ** 2)\n",
    "    y_mean = np.mean(y_neigh)\n",
    "    ss_tot = np.sum((y_neigh - y_mean) ** 2)\n",
    "    \n",
    "    # 处理分母为零的情况\n",
    "    if ss_tot == 0:\n",
    "        r2 = 0.0\n",
    "    else:\n",
    "        r2 = 1 - (ss_res / ss_tot)\n",
    "    local_R2.append(r2)\n",
    "\n",
    "\n",
    "# 将结果保存到DataFrame\n",
    "output_data = data[['FID_','ENG_NAME','NAME_1','X_position', 'Y_position',\n",
    "                    'SMCI_2020', 'GDP_2020', '人口_2020','土地利用结构_2020'\n",
    "          ,'人口密度_2020','年均温_2020','人均GDP_2020','年降水量_2020',\n",
    "          '燃烧面积占比_2020_2','火点个数_2020','单位面积FRP_2020']].copy()\n",
    "\n",
    "# 输出残差、预测值、local_R2\n",
    "output_data['residual'] = results.resid_response\n",
    "output_data['y_hat'] = results.predy\n",
    "output_data['local_R2'] = local_R2\n",
    "\n",
    "\n",
    "# 添加各变量的系数和t值\n",
    "for i, col in enumerate(['截距','SMCI_2020', 'GDP_2020', '人口_2020','土地利用结构_2020'\n",
    "          ,'人口密度_2020','年均温_2020','人均GDP_2020','年降水量_2020',\n",
    "          '燃烧面积占比_2020_2','火点个数_2020','单位面积FRP_2020']):\n",
    "    output_data[f'{col}_coefficient'] = results.params[:, i]\n",
    "    output_data[f'{col}_t'] = results.tvalues[:, i]\n",
    "\n",
    "# 保存为CSV\n",
    "output_data.to_csv('D:\\Lenovo\\Desktop\\云南大学\\毕业设计\\毕设数据\\输出数据\\csv\\MGWR结果\\MGWR_2020_3.csv', index=False,encoding='GB2312')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shapley分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(72,)\n",
      "Provided model function fails when applied to the provided data set.\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "'numpy.ndarray' object is not callable",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[17], line 23\u001b[0m\n\u001b[0;32m     20\u001b[0m \u001b[38;5;28mprint\u001b[39m(y_test_scaled\u001b[38;5;241m.\u001b[39mshape)\n\u001b[0;32m     22\u001b[0m \u001b[38;5;66;03m# 5. 计算SHAP值\u001b[39;00m\n\u001b[1;32m---> 23\u001b[0m explainer \u001b[38;5;241m=\u001b[39m \u001b[43mshap\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mKernelExplainer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmodel\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43my_test_scaled\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mX_train_scaled\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     24\u001b[0m shap_values \u001b[38;5;241m=\u001b[39m explainer\u001b[38;5;241m.\u001b[39mshap_values(X_test_scaled,gc_collect\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[0;32m     26\u001b[0m \u001b[38;5;66;03m# 6. 绘制beeswarm图\u001b[39;00m\n",
      "File \u001b[1;32m~\\AppData\\Roaming\\Python\\Python311\\site-packages\\shap\\explainers\\_kernel.py:97\u001b[0m, in \u001b[0;36mKernelExplainer.__init__\u001b[1;34m(self, model, data, feature_names, link, **kwargs)\u001b[0m\n\u001b[0;32m     95\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmodel \u001b[38;5;241m=\u001b[39m convert_to_model(model, keep_index\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mkeep_index)\n\u001b[0;32m     96\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m convert_to_data(data, keep_index\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mkeep_index)\n\u001b[1;32m---> 97\u001b[0m model_null \u001b[38;5;241m=\u001b[39m \u001b[43mmatch_model_to_data\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmodel\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     99\u001b[0m \u001b[38;5;66;03m# enforce our current input type limitations\u001b[39;00m\n\u001b[0;32m    100\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata, (DenseData, SparseData)):\n",
      "File \u001b[1;32m~\\AppData\\Roaming\\Python\\Python311\\site-packages\\shap\\utils\\_legacy.py:142\u001b[0m, in \u001b[0;36mmatch_model_to_data\u001b[1;34m(model, data)\u001b[0m\n\u001b[0;32m    140\u001b[0m         out_val \u001b[38;5;241m=\u001b[39m model\u001b[38;5;241m.\u001b[39mf(data\u001b[38;5;241m.\u001b[39mconvert_to_df())\n\u001b[0;32m    141\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 142\u001b[0m         out_val \u001b[38;5;241m=\u001b[39m \u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    143\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m:\n\u001b[0;32m    144\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mProvided model function fails when applied to the provided data set.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[1;31mTypeError\u001b[0m: 'numpy.ndarray' object is not callable"
     ]
    }
   ],
   "source": [
    "# 分割训练集和测试集数据\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.2, random_state=42\n",
    ")\n",
    "def predict(FID_):\n",
    "\n",
    "    return y_hat\n",
    "\n",
    "# 计算SHAP值\n",
    "explainer = shap.KernelExplainer(model=predict(X_train[\"FID_\"]), data=X_train.drop[\"FID_\"])\n",
    "shap_values = explainer.shap_values(X_test,gc_collect=True)\n",
    "\n",
    "# 绘制beeswarm图\n",
    "plt.figure(figsize=(10, 6))\n",
    "shap.plots.beeswarm(shap_values, show=False)\n",
    "plt.title(\"SHAP Value Beeswarm Plot\")\n",
    "plt.xlim(-1,1)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下是测试内容"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(results.adj_alpha)\n",
    "# print(results.influ is None)\n",
    "# print(results.adj_alpha_j)\n",
    "# print(results.resid_response.shape)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gdal-env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
